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The Anderson model for a single impurity coupled to two leads is studied using the 
GW approximation in the strong electron-electron interaction regime as a function 
of the alignment of the impurity level relative to the chemical potentials in the 
leads. We employ a non-equilibrium Green's function technique to calculate the 
electron self-energy, the spin density and the current as a function of bias across the 
junction. In addition we develop an expression for the change in the expectation 
value of the energy of the system that results when the impurity is coupled to the 
leads, including the role of Coulomb interactions through the electron self energy 
in the region of the junction. The current-voltage characteristics calculated within 
the GW approximation exhibit Coulomb blockade. Depending on the gate voltage 
and applied bias, we find that there can be more than one steady-state solution for 
the system, which may give rise to a hysteresis in the I-V characteristics. We show 
that the hysteresis is an artifact of the GW approximation and would not survive if 
quantum fluctuations beyond the GW approximation are included. 



I. INTRODUCTION 



Transport through nanoscale junctions poses a number of interesting physical problems. 



In particular, electron-electron interaction effects may be important, as evidenced by the 
observation of phenomena such as the Coulomb blockade and the Kondo effect [l, 3]. The 
local electronic structure is also important. The energy and character of the electronic states 
in the junction region that are responsible for electron transport will depend on the details 
of bonding between the molecule and the electrode. This has motivated the use of ab initio 
theories for electron transport through nanostructures that are based on Density Functional 
Theory (D 
properly [3. 

the level of model systems, a complete solution of the nonequilibrium interacting electron 



T). However, local density functionals do not treat the discreteness of charge 
J]. In particular Coulomb blockade phenomena become problematic. Even on 



problem is not available. The numerical methods which work so well in equi 



beginning to be applied to non-equilibrium systems 
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ibrium are only 
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oring selfconsistent perturbative and other, nonperturbative approaches [13 
3, Q, Q However, a complete treatment which can be extended to incorporate 



actual, junction-specific aspects is not yet available. 

In this work, we study a model system, namely the single impurity Anderson model 



2l| coupled to two leads. We use a Green's function approach to calculate the properties 



of the junction, both in equilibrium and as a function of applied bias across the junction. 
The electron-electron interactions are incorporated through the electron self energy operator 
on the impurity, using an out-of-equilibrium generalization of the GW approximation 22]. 
Using this approach we can calculate the local spin density in the junction and the current 
as a function of bias. In addition we develop and apply an extension to non-zero bias of the 
usual expression 23) for the change in the average energy of the impurity due to coupling to 



the leads. The GW approximation has been widely and successfully used to study e 



excitations in materials at equilibrium with a realistic, atomic scale description 24 



ectronic 
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28]. This is one of the motivations to study the out-of-equilibrium generalization for 

In particular, the intermediate coupling/interaction 



nanoscale junctions 



13, Q 



regime of the sing le impurity Anderson model has recently been studied using the GW 
approximation 

We are interested in the intermediate to strong coupling regime, in which Coulomb block- 



ade effects are important. At equilibrium and for zero temperature, as the local Coulomb 
interaction on the impurity is increased (relative to the hybridization with the leads) a lo- 
cal moment forms. In the limit of ksT — ► and vanishing bias, the local moment on the 
impurity is quenched through formation of a singlet ground state. The spectral function 
splits into three parts, two Hubbard bands and one central Kondo peak. In a closely re- 
lated earlier study [29], it was shown that in the regime of intermediate strength of the 
Coulomb interaction, the GW approximation provides an incorrect representation of the 
linear response conductance. In fact, this regime is not well described at equilibrium even 



by more sop 



imation 



listicated perturbative approaches, such as the fluctuation-exchange approx- 



30, 



311 ] . Here we probe the strong coupling, Coulomb blockade regime. In this 



regime, the Kondo temperature Tk becomes very small and at experimentally relevant tem- 
perature scales the Kondo peak will be washed out. Similarly, when considering bias large 



compared to the Kondo temperature, the Kondo peak also gets washed out [32|, |33j. In 
these regimes a self-consistent perturbative approach may be adequate. We find through 
non-equilibrium calculations that the self-consistent GW approximation can describe impor- 
tant features of the Coulomb blockade regime, such as the Coulomb diamond signature with 



ecule transistors 

3- 



no Kondo-assisted tunneling, in accordance with experiments on single-mo 
characterized by weak effective coupling between molecule and electrodes 

The non-equilibrium GW calculations exhibit hysteresis in the IV characteristics: at 
some values of applied bias and gate voltage, there is more than one steady state solution. 
A related example of bistability has been found in DFT calculations of a junction involving 
an organometallic molecule 3_4j. However, we believe that in the problem that we study 
here, the hysteresis is an artifact of the approximation [3] . In fundamental terms, a molec- 
ular junction is a quantum field theory in space and 1 time dimension. Model system 
calculations 



18 



191 ] have confirmed that departures from equilibrium act as an effective 



temperature which allows the system to explore all of its phase space, preventing bistability 
from occurring. We will show by an energy calculation that in the present problem similar 
processes exist. 

The rest of the paper is organized as follows. In Section II, the model Hamiltonian is de- 
scribed. Section III presents the non-equilibrium, self-consistent Green's function approach 
that we use, including the GW approximation, an expression for the change in the average 
energy as well as an expression for the current that allows to distinguish the Landauer-like 



and the non-coherent contributions. The results of the calculations for the single impurity 
Anderson model are developed in Section IV. Derivations of the expressions for the physical 
observables appear in Appendices A, B and C. 

II. MODEL HAMILTONIAN 

We consider the Anderson model for an impurity coupled symmetrically to non- 
interacting leads. We are interested in steady-state solutions of this system. The Hamil- 
tonian describing the system, H, can be written as a sum of a non- interacting part, Hq, 
plus an interacting one, H e _ e , describing the electron-electron interaction in the impurity: 
H = H + H e _ e . 

The non-interacting part is treated at the tight-binding level (Fig. The left (L) 

and right (R) leads are modeled as semi-infinite chains of atoms (i=l,... oo or -1,... — oo), 
characterized by the hopping parameter t and chemical potentials /i^ and /iR. We choose 
t = 5, resulting in the band-width of the metallic leads extending to ±10 about the chemical 
potential of each lead which we fix at the center of each electrode band. The system is driven 
out of equilibrium by applying a source-drain bias voltage V, setting fi L = —fiR = V/2; 
the impurity levels can also be shifted according to a gate voltage Vq (Fig. db). The 
hybridization term describes the coupling between the impurity (site 0) and the nearest 
atoms of the two leads (sites ±1), and is parameterized according to the hoping parameter 

7- 

-2 oo 

H = u L N l +h r N r +V g n -t( + E) 5^(cl c i+ia+4+i CT c ^)-7 E ^2( c l c 0a+cl a c ia ) 
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(1) 

where Num are the electron number operators in the L(K) leads: 

— 1 oo 

n l = E c L c - ^ = E E 4*- (2) 

i=— OO (7 i=l IT 

and no is the electron number in the impurity: 

n o = E c o° c °° ( 3 ) 

a 

The electron-electron interaction inside the impurity is taken into account through the 



usual U-term: 



U %t™0| = \ Yl C 0, a 4,/3^«',/3/3' C 0,/3'C0,a' (4) 



H e - e . „,. 2 

a,a',p,p 



There are several choices we can make for the 2-particle interaction V aa ',pp'- We choose one 
that describes non-spin-flip scattering: 



and has a spin-dependent form: 



V aa ',pp> — V a @ 5 aa i 5pp> (5) 



V a p = U (1 - 5 a p) (6) 



Another choice for the 2-particle interaction, which results in the same Hamiltonian as in 
Eq. (pEJ), would be one with a spin- independent form: V a p = U . However, in the context 
of the GW approximation for the Anderson model, the spin-dependent form is a better 
choice (29!. Indeed, it has been shown that the spurious self-interactions can be a major 
source of error in transport calculations, especially when the coupling to the leads is weak 
4j. Comparing the two choices for V a p, the spin-dependent one has the advantage of being 
free of self-interaction effects, and it also accounts for more quantum fluctuations in the 



spin-spin channel 



29]. 



In the present model, the potential due to the applied source-drain bias V and gate 
voltage Vq changes only at the junction contacts (Fig. [Ib). Also, the direct electron- 
electron interaction between the impurity and the leads is neglected. These approximations 
are justified in realistic systems in which the screening length in the leads is very short. 
We shall be interested in the limit of very small effective coupling to the leads T = jt. 
Our choices 7 = 0.35 and t — 5 imply T = 0.05. The on-site Coulomb repulsion between a 
spin-up and a spin-down impurity electron is set to U — 4.78 ~ 100r. At equilibrium and 
half-filling, the Kondo temperature Tk is then [35 ] : 



Tk ~ 0.2v / 2fZ7exp(-7rf//8r) (7) 



which is thus negligible small. The results that we present for this set of parameters hold, 
qualitatively, for a wide range of parameters consistent with a weak hybridization and strong 
Coulomb interaction regime. 



III. SELF-CONSISTENT NON-EQUILIBRIUM GREEN'S FUNCTION 

FORMALISM 



A. Hamiltonian and Basic Formalism 



Electron correlation effects in the impurity are studied using a non-equilibrium Green's 
function formalism, by solving self-consistently for the variou s [r etarded (r), advanced (a), 
lesser (<) and greater (>)] Green's functions of the impurity j^g, 



G r (u) = [(u - V G )I - A£M - A^cu) - V H - TI\u)}- 1 
G<(cu) = G r {u)[if L {u)Y L {u) + if R {u)Y R {u) + £<H]G a H 



(8) 
(9) 



where all quantities are matrices in the space spanned by the junction degrees of freedom, 
in the present case the up and down components of the impurity spin 381 ] . 

Above, A' r stands for the retarded lead self-energy, which, for our model Hamiltonian, 
takes the form 391 ] : 
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(10) 



and we have used the notation: 



Yl( R )(uj) = i' Ar 



i[A 



L(R) 



V)-A 



L(R) 



UJ 



,tl 



(11) 



The hybridization functions Al(r) are centered on the chemical potentials ^l(r), such that 
the isolated leads are neutral. 

V represents the Hartree potential: 



—II J 



(12) 



and Y7 (S < ) is the retarded (lesser) impurity self-energy, describing the effects of electron- 
correlation inside the junction. The electron occupation numbers appearing in Eq. ([9]) are 
the usual statistical factors for a system of electrons: Jur) (u) = l/{exp[(u — Hl{r))I ksT] + 
1}. Since we operate in the regime of very small Kondo temperature, we envision choosing 



an experimentally relevant temperature that is large compared to T%, but which is much 
smaller than the coupling to the electrodes. 

The other two non-equilibrium impurity Green's functions can be simply obtained using: 



G>(u) = G r (uo) - G a (uo) + G < (uo) 



(13) 
(14) 



B. The GW approximation for the impurity self-energy 

In the GW approximation for the electron self-energy, one does perturbation theory in 
terms of the screened interaction W, keeping the first term in the expansion, the so called 



GW diagram. The GW approximation has long been successful 



equilibrium quasiparticle properties of real materials 24 
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y used in describing the 
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28]. It has also been 



applied to the study of real materials out of equilibrium, such as highly irradiated semicon- 
ductors J40I] . or, more recently, in transport calculations through molecular nanojunctions 
[14] . [lol . 17]. For equilibrium properties the GW approximation has been compared to a 
numerically exact quantum Monte Carlo treatment 29J]; it has been found to be adequate 
for small interactions or for high T, but not in the mixed valence or Kondo regimes. 

Within the out-of-equilibrium GW approximation, the general self-energy expressions 



have the following form in frequency space 



40|: 



KA") 



2^ g<ae)w;A" - e) + ij — g:ae)w>A" - m ^) 

=ij^ G< a ,(E)W<,(u - E) (16) 

where the screened interaction W can be obtained from the irreducible polarizability P 
through: 

W r {uj) = [I- VP^u)}-^ (17) 
W < {u) = W r {u)P < {u)W a {u) (18) 
W>(u) = W r (oo)P > (uj)W a (uj) (19) 
The irreducible polarization P is evaluated in the random phase approximation (RPA): 

JTp 

J G< a ,{E)Gl, c {E-u) (20) 
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Jt^ G «AE) G< a (E 



2tt 



— GIAE) G< a (E - u )-iJ — G< a ,(E) G r a , a (E - u) (21) 
% J ^ G «AE) G> a (E - uj) (22) 



Setting P = yields the Hartree-Fock approximation. 

The set of equations for G, S, W and P are solved to self-consistency, starting from an 
initial condition for G. All the quantities are calculated on a real frequency grid (either 
regular or log-scale), with an cu-range up to ±10t. Real and imaginary parts of the vari- 
ous quantities are calculated explicitly, making sure that the retarded functions obey the 
Kramers-Kronig relation. In order to speed up the self-consistent process, we employ the 
Pulay scheme to mix the Green's functions using previous 

iterate solutions 
Gt: 1 = (1 - a)Ql + aQi ut (23) 
where Q n are constructed from the previous m iterations: 

m 
i=l 

and we choose three components for the parameter vector Q: 9ftG r , ^sG r and < ^sG < . The 
values of $ are obtained by minimizing the distance between Q\ n and Q J out - The scalar 
product in the parameter space is defined using the integral in Fourier space of a product 
of the component Green's functions. We found the speed of the convergence process to be 
quite independent on the choice of reasonable values for m, as well as on the number of 
components for the parameter vector Q. As for the parameter a, smaller values (< 0.1) were 
needed for small bias voltages (V < 0.5), while a = 0.4 was sufficient in order to achieve 
fast convergence for larger biases. 

C. Relation to physical observables 

The Green's functions of the impurity can be used to extract information about observ- 
ables pertaining to the impurity or even to the leads. Thus, the spectral function of the 
impurity A(u) is simply related to the retarded Green's function: 

A(uj) = — Tr%G r (uj) (25) 

7T 



where Tr stands for trace over the impurity spin degrees of freedom. Also, the average 
impurity spin occupation number is: 

M = J ^ G<» (26) 

The expression for the average current passing through the junction is given by the general 
Meir-Wingreen expression 42j, which can be recast as (see Appendix A for the derivation): 



/ = J dw [f L (u) - f R (u)} Tr{T L (u) G r (cu) T R (cu) G»} 

+ Jdu Tr{[T L (u)-r R (u)} G r (u) [±E<(u>)] G»} 

+ J duo Tr{[f L {u)Y L {u) - f R (uj)T R (uj)} G» G a (oo)} (27) 

The first (Landauer type) term plays an important role whenever correlations beyond the 
Hartree-Fock level are not considerable. It gives the coherent component of the current. The 
second term is in general very small for symmetric leads with relatively wide bands, when 
Tl(uj) ~ T r (uj). The last term becomes important when the electron-electron correlation 
effects are such that — 3£ r « ^l{r)- 

Having an expression for the average energy associated with the junction for non- 
equilibrium can be useful for a number of purposes, including calculation of current de- 
pendent forces 43J]. By formulating this as the difference 5S between the average energy of 
the total system (leads coupled to impurity) and the average energy of the isolated leads, 
a finite result can be obtained. This can be done starting with the following expression for 
the total average energy of the system [46 ]: 

£ = \ f^Tr{(H + uI)G<(u)} (28) 
2 J Am 

where the trace Tr is taken over a complete set of states spanning the junction (indices n) 
and the leads (indices k). Alternatively, an equation of motion approach can be used 23]. 
We find that the two approaches give the same results. The first approach is presented in 
Appendix B. Naturally, the energy can be decomposed into three terms, related respectively 
to the average energy of the impurity ^ mp , the average energy of interaction between leads 
and impurity £im P -i e ads, and the average energy difference in the leads before and after 
adding the impurity 5£ leads : 

£imp ")~ £imp— leads $ Pleads (^) 



where: 



with: 



'-'imp 



duj 
2tH 



(lu + Vg) TrG<{uo) 
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imp— leads 



duj 



Tr{^l{u>)+MSr R {u>)] G<{u) 



(30) 



(31) 



68, 



leads 



- i[f L {u) ZAKtu) + f R {u) 9fAi(w)] [G» + G»]} 
p:Tr{&F L {u>) + toF R {u>)] G<(co) 

2,711 

- i[ f L (u) ZF L (cu) + f R (u) %F R (u;)} [G a (uj) + £»]} (32) 



2oj-^A r T(R) 
duj L{R) 



F L (R )nm (w) = - A r L(R) (w) - 2^— Ai (J0 (33) 

We note that the average energy change in the two leads is always finite in the steady state 
case. A similar statement holds for the average number of electrons displaced in the two 
leads 6Ni eads (explicit expression in Appendix C). 



IV. RESULTS 



A. Coulomb blockade 



In the weak coupling/strong interaction regime, the electron transport through a junction 
can be blocked due to the charging energy in the junction. Figure [2](a) shows the calculated 
impurity occupation number (no) = (no\) + ( n oi) as a function of the gate voltage Vg, at 
zero applied bias V = 44J. One can clearly see the Coulomb staircase. The electron-hole 
symmetry of the Hamiltonian describing the system, H, insures that the spectral function 
satisfies: A(u; Vg+U/2) — A(— ui; — Vq — U/2). As a consequence, one has: {ti (Vg + U/2)) = 
2 — (n (— Vg — U/2)). A similar Coulomb staircase picture can be obtained at the Hartree- 
Fock approximation level. 

The impurity occupation number evolves from to 2 as Vg is decreased from positive to 
negative values. Figure Efb) shows the evolution of the spectral function for three repre- 
sentative values of Vq- For Vg + U/2 = ±4, the solution is non-magnetic, with both spin 
levels degenerate, empty or occupied. At the symmetric point (half- filling) Vq+U/2—0, the 
solution is a broken symmetry magnetic ground-state, with one spin occupied and the other 
empty. Since we consider temperatures that, although small, are still large compared to Tk, 



the degenerate magnetic ground state is an appropriate representation of the physics. In 
Fig. [2(a), the magnetic solution is found for \V G + U/2\ < 2; for 2 < \V G + U/2\ < 3, a well 
converged (non-magnetic) solution could not be found at k^T = 0. 

Figure [3(a) shows a color-scale plot of the current I as a function of the applied bias V 
and gate voltage V G . The plot is obtained by forward scan of the bias, i.e. using the lower 
bias solution as starting input for the higher bias calculation. One can see the formation 
of Coulomb diamonds, inside which the current is negligible, a signature of the Coulomb 
blockade regime. A similar color-scale plot of the differential conductivity would show sharp 
peaks at the edges of the Coulomb diamonds, but no tunneling channel in the zero bias region 
inside the central Coulomb diamond. Such a tunneling channel is absent in experiments on 
single-molecule transistors characterized by weak coupling between molecules and electrodes 
but has been observed when coupling to the electrodes is strong enough that the Kondo 
temperature is appreciable T K ~ 10 to 30 meV [l,^]. 

At zero bias, zero temperature and at the symmetric point, the unitarity limit 48] requires 
that the differential conductivity equals 2e 2 /h. The broken (magnetic) symmetry solution 
in the GW approximation in the strong interaction regime does not satisfy the unitarity 
limit; the spectral function does not have the correct height near the chemical potential. 
Therefore, the GW approximation can not account for the zero-bias tunneling channel ob- 
served for T < T%. Under finite bias, the differential conductance due to the Kondo peaks 
in the spectral function must fall off once the bias exceeds the Kondo temperature 32|; the 
Kondo peak splits under non-zero bias, following the two different chemical potentials and 
broadens quickly with increasing bias. Therefore, the width in applied bias for which such 
a channel would be observed in the exact theory is of order Tk- For the strong interaction 
regime considered here, this is negligible. Thus, the GW approximation provides the correct 
qualitative features of the Coulomb blockade regime, namely Coulomb diamonds with no 
Kondo-assisted conductance channels. 

The size of the Coulomb diamond depends on the interplay between the repulsion U and 
the coupling to the leads, T. In the limit of U/Y — > oo the system becomes effectively 
an isolated ion, and the size of the diamond is set by U. In our case, U/Y « 100 and 
the computed size of the Coulomb diamond is only slightly smaller (by ~ 20%) in the 
GW approximation than in the Hartree-Fock approximation. However, we suspect that the 
magnetic solution found in the GW approximation underestimates the electronic correlation 



originating from spin-spin quantum fluctuations, and thus a more exact theory should result 
in smaller size Coulomb diamonds than the ones we find. 

The corresponding average electron occupation number (no) is shown in Fig. [3](b), where 
we can see that (no) takes integer values of 0, 1 and 2 inside the Coulomb diamonds. For 
a given gate voltage, the spectral function of the system changes appreciably only when 
the left or right lead Fermi levels get closer to one of the impurity resonance levels. As 
soon as a resonant level is pinned by a Fermi level, the current increases while the impurity 
occupation number either increases or decreases depending whether the pinned level is empty 
or occupied. 



B. Hysteresis in the I-V characteristics 



In an earlier study 22| we concluded that, in the regime of intermediate strength of the 
Coulomb interaction, the GW approximation leads to a broken spin symmetry ground state 
and thus fails to describe the spectral function correctly, missing completely the Kondo 
peak. A non-magnetic solution in the interaction regime U/T > 8 has been elusive for 



other authors as well 16|. Recently, by employing a logarithmic frequency scale near the 
Fermi level, we have been able to find a non-magnetic solution in the strong interaction 
regime up to U/T « 25 and ksT = 0. Our results [45] show that equilibrium properties 
of the Anderson model, such as the total energy, Kondo temperature, T-linear coefficient 
of the specific heat or linear response conductance, are not satisfactorily described by the 
non-magnetic solution in the GW approximation, as it was previously noted for several of 
these properties 0, 

For the interaction strength considered in the present work, U/T ~ 100, we have been 
able to calculate the non-magnetic solution at zero bias by considering small non-zero tem- 
peratures. We will consider fc#T = 0.01 throughout the rest of the paper. Figure BJa) 
shows the impurity occupation number as a function of gate voltage for the non-magnetic 
solution. We see that the Coulomb blockade plateau is not properly described; the impurity 
occupation number changes linearly about the symmetric point Vq + U/2 — 0. Figure H^b) 
shows the spectral function associated with the non-magnetic solution for two representative 
cases. In the symmetric case, one sees a broad peak (whose width is set by U) with a nar- 
row portion near E = (whose width is set by ksT). As the gate voltage is changed from 



the symmetric point, the narrow portion remains pinned near E = 0, but the broad peak 
shifts together with Vq, hence the linear change in < n occ > as observed in Fig. H](a). Near 
Vg + U/2 = ±2.8, the non-magnetic solution cannot sustain a narrow portion near E — 0, 
and the solution jumps into a phase with one narrow peak (of width ~ V) away from E = 
(as seen in Fig. E](b) for Vq + U/2 — 4). In the region of gate bias near the transition at 
Vq + U/2 = ±2.8, the calculations get increasingly difficult to converge; for some values of 
the gate voltage a converged solution with retarded functions obeying the Kramers-Kronig 
relation could not be found. 

Figure EJ^a) shows the current through the junction / as a function of the applied bias V 
for a specific gate voltage Vg, Vg + U/2 — such that the system is at half-filling, (n ) = 1. 
The general results do not depend on this symmetry. The same qualitative results hold for a 
broad range of gate bias \ Vg + U/2\ < 2. Results obtained both in the GW and Hartree-Fock 
approximation are shown. These include a forward scan, starting from zero applied bias, and 
a reverse scan starting from V = 8. At zero bias, we start with the magnetic solution, with 
one spin level occupied and the other one empty, as shown by the solid-line curve of Fig. 
Mb). Then in the forward scan, the initial input at higher bias is taken from the converged 
solution at lower bias. For the reverse scan, the opposite approach is taken. Note that the 
use of ksT = 0.01 has essentially no effect on the results except for the reverse scan with 

V < 0.5 where the finite temperature helps to stabilize the self consistent magnetic solution. 
Also, for reference, the I — V data shown in Fig. E^a) was obtained by forward bias scan. 

In both the Hartree-Fock and the GW approximation, as the bias is increased, the two 
spin levels remain outside the bias window and the current is negligible until V approaches 
a value of order (but less than) U. In Hartree-Fock this value is V ~ 4.0, while in GW it is 

V ~ 3.2. At this point, where the broadened impurity levels get pinned by the two chemical 
potentials, the character of the steady-state solution changes from magnetic to non-magnetic. 
At this bias, the current increases suddenly. Correspondingly, the spectral function shows 
one double-degenerate peak centered half-way in between the two chemical potentials (Fig. 
EJ^b)). For higher bias, the Hartree-Fock and GW approximations result in qualitatively 
and quantitatively different behavior. In Hartree-Fock, the current is approximately pinned 
at the value expected for a single, half-filled resonance in the bias window (2irTe/h). The 
overall downward drop is explained by the finite band width of the electrodes. However, in 
the GW approximation, the spectral function shows substantially larger broadening and the 



current increases steadily with bias as the spectral weight inside the bias window increases. 
Correspondingly, upon analysis of contributions to the current in this regime, it is largely due 
to non-coherent transport, as — 3S r >> T and the main component of the current is given by 
the last term of the right hand side of Eq. ( 1271) . The backward bias scan is started from the 
non-magnetic solution at V — 8. As the bias is decreased, the solution remains non-magnetic 
well below the transition bias point from the forward bias scan, resulting in hysteresis in the 
I — V curve. While in Hartree-Fock, the current remains high down to relatively low applied 
bias, the calculated current in the GW approximation drops approximately linearly. 

The physical description of the magnetic solution is straightforward. The spectral func- 
tion shows two peaks (Fig. E](b)), spin up and spin down, one occupied and the other one 
empty, separated in frequency by a little less than U. The results from the GW approxima- 
tion are very close to those from Hartree-Fock approximation in this case. There are very 
few occupied-to-empty electron-hole same-spin excitations; the polarization P is very small. 

The non-magnetic solution is more complex and the physical picture is rather different 
for the Hartree-Fock and the GW approximations. While for Hartree-Fock the spectral 
function showing only one sharp peak with width equal to T, the spectral function in the GW 
approximation is much broader (Fig. [5](b)). While the overall broadening depends strongly 
on the interaction parameter U, the applied bias V affects the region of width V about 
E — (Fig. E(c)). Furthermore, the width of the spectral function is almost independent 
of the effective coupling coefficient V. For example, for ksT = 0.01 and V G [0,8], the 
spectral function plot for T = 0.1 is almost undistinguishable from that for the T = 0.05 
case. This indicates that the broadening is due to quantum fluctuations taking place on the 
impurity. The applied bias dependent broadening can be traced back to the large imaginary 
part of the retarded self-energy, as shown in Fig. El At zero bias and zero temperature the 
Fermi liquid behavior of the system guarantees 3S r (0) = 0. The non-zero value of ImS r (0) 
shown in Fig. [6] is clearly a non-equilibrium, non-zero bias effect. A similar broadening, 
increasing strongly with bias, has been also observed in recent calculations based on the 
GW approximation for a two- level model molecule 171 ]. 

The broadening of the spectral function for the non-magnetic solution in the GW ap- 
proximation can be understood by looking at how the spectral function and the retarded 
self-energy changes as we iterate the non-magnetic solution from Hartree-Fock to GW. Here 
we denote with GqWq the intermediate solution obtained with the Hartree-Fock Green's func- 



tions as input. At the Hartree-Fock level, the non-magnetic solution has one narrow central 
peak, with half-width at half- maximum approximately given by T = 0.05. The entire peak 
is situated inside the bias window, as shown in Fig. [3(a). In that energy range one can find 
both occupied and empty (more exactly half-occupied) quasi-states of the same spin. Now, 
such a quasi-state can easily decay into another quasi-state with lower or higher energy, 
by emitting or absorbing an electron-hole same-spin excitation with energy within the bias 
window range. Thus, 3X7, which is proportional to the inverse lifetime of the quasi-state, 
becomes very large at the the GqWq level, as seen in Fig. [3(b). From the G W result for 
3?S r (related to SE'' through a Kramers-Kronig relation), it follows that the GqWq spectral 
function shows two double degenerate (spin up and spin down) peaks, situated outside the 
bias window (Fig. [7(a)). If, at each of the next iterative steps i, we would use as input only 
the Green's functions from iteration i — 1, the spectral function would oscillate between the 
two types (Hartree-Fock and GqWq) of solution. However, by means of the Pulay mixing 
scheme, we are able to achieve convergence rather fast, with the self-consistent GW solution 
looking somehow in between Hartree-Fock and GqWq, as seen in Fig. [5(b). 

The calculated I-V curves in Fig. [5](a) result from the existence of two steady state 
solutions over a broad range of applied bias that are accessed depending on initial conditions. 
Our procedure of stepping the applied bias in forward followed by reverse scans with self 
consistent solution at each step simulates an adiabatic voltage scan and the existence of two 
stable solutions results in hysteresis. One may ask whether quantum fluctuations that are 
beyond the scope of the GW approximation would eliminate the hysteresis. To probe this, 
we need to understand the energy difference between the system in the magnetic and the 
non-magnetic solutions in the hysteretic region. Figure [8] shows the change in the average 
energy of the total system, 58, calculated as described in Section III(B), as a function of 
the applied bias at half-filling. Results are shown for both the Hartree-Fock and the GW 
approximations, following the same loop of forward and reverse bias scans. For weak effective 
coupling between impurity and leads, for the magnetic solution, one has SS ~ £i mp + 0(T). 
Near equilibrium, the magnetic solutions in the forward bias scan show very similar energies, 
close to the energy of the isolated, single-occupied impurity: 5£ mag ~ Vq — —U /2. However, 
at the applied bias where the current rapidly increases and the solution changes to non- 
magnetic, Hartree-Fock yields an average energy higher than the magnetic one by about 
{7/4. On the other hand, the GW approximation shows an average energy change that is 



much smaller. Correspondingly, on the reverse bias scan, the bias dependence of the average 
energy is also much different. While the energy in the Hartree-Fock approximation remains 
high as the bias approaches zero, the energy in the GW approximation approaches a value 
that is only higher than the zero bias magnetic state by about T/20. 

We have found that in the strong interaction regime, there are two distinct self consistent 
solutions with the GW approximation. These lead to hysteresis in the calculated I — V 
curves. However, at zero bias, bistability is forbidden for the Anderson model Q]- 
Therefore, the states represented by those solutions found in the GW approximation must 
be unstable with respect to quantum fluctuations that have not been taken into account. 
The fact that the average energy of the magnetic state is lower than that of non-magnetic 
solution, is probably an indication of the larger weight of the magnetic solution in the 
emerging exact many-body state. As the bias is increased away from equilibrium, Fig. [H] 
shows that the energy difference between non-magnetic and magnetic configurations also 
increases in the GW approximation. However, for applied bias larger than about T/20, the 
energy difference is smaller than the applied bias. This means that at non-zero biases on-shell 
processes will be possible through which one configuration can decay into the other one (with 
one electron transferring from one lead to the other to insure total energy conservation). We 
thus expect that out-of-equilibrium, the lifetime of the GW bistable states would be even 
smaller than at equilibrium. Quantum fluctuations between the two degenerate magnetic 
configurations and the non-magnetic one will eliminate the hysteresis and renormalize in a 
non-trivial way the emerging unique many-body state. Therefore, the hysteresis in the I — V 
curve is probably another signal that the GW approximation is not representing important 
aspects of the strong interaction regime. A calculation of the lifetime of the bistable states 
found with the GW approximation is beyond the scope of the present work, but would be 
very valuable. 



V. SUMMARY 



In this work we used the GW approximation to study the role of electron-electron corre- 
lation effects in the out-of-equilibrium single impurity Anderson model. We considered the 
regime with weak level broadening and strong Coulomb interaction, treating the electron- 
electron interaction with the self-consistent GW approximation for the electron self energy. 



We found that the GW approximation accounts for Coulomb blockade effects. The low 
conductance (blockade) region in gate bias and source-drain bias corresponds to a magnetic 
solution in the GW approximation. At the edge of the blockade region, the current jumps 
and the self consistent solution changes to a non-magnetic character. The position of the 
transition and the jump in current are renormalized from the Hartree-Fock values. However, 
we also found a self consistent non-magnetic solution inside the Coulomb blockade region. 
As a consequence, the GW approximation also predicts an unphysical hysteresis in the I- 
V characteristics of the system. Outside the blockade region, e.g. where the source-drain 
bias is high and the magnetic solution is not stable, we expect that the GW approximation 
gives a reasonable account of the conductance. However, the jump in current at the edge of 
the blockade region and the hysteresis inside the blockade region both appear to arise from 
a first-order-transition-like bistability in the GW approximation. An analysis of the total 
energy difference between the magnetic and non-magnetic solutions suggests that quantum 
fluctuations beyond the scope of the GW approximation would result in rapid decay of the 
non- magnetic solution, eliminating both the sharp jump and the hysteresis. 
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APPENDIX A: AVERAGE CURRENT THROUGH THE JUNCTION 



We start from the Meir-Wingreen expression for the current from the left lead (in units 



of e = h= 1) [42|: 

h = i J du Tr{T L (u) G<(u) + f L (w)r L (u) [G r (u) - G a (u)}} = J duJ L {u). (A-l) 
Using that in steady state I = II = (II — /r)/2, and making use of Eq. (jHJ) and the relation: 
G r (u) - G a {u) = G r (cu) [A r L (u) + A r R (cu) + S r (cu) - h.c] G a {u) (A-2) 

one obtains: 

1 = 1 [du [f L (u) - f R (u)] Tr{Y L {u) G» T r (lu) G»} 

+ \j du [f L (co) - f R (u)] Tr{T R (u) G r {u) T l (lu) G»} 
+ \J du Tr{[T L (uj) - r R (u)} G r (u) E<(w) G a (co)} 

+ ~ J doo Tr{[f L (u)T L (u) - f R (u)r R (u)] G r (u) [E» - £»t] G»} (A-3) 

In the single impurity Anderson model case, the Green's functions are symmetric (the off- 
diagonal elements being simply zero) and the first two terms in Eq. f]A-3|) are equal, with 
the final expression for the current reading as in Eq. (127)) . 

APPENDIX B: THE CHANGE IN ENERGY CAUSED BY IMPURITY 

For simplicity, we consider eigenstates of the non-interacting isolated junction (energies 
e n ) and isolated leads (energies Denoting with g the Green's function of the isolated 
lead, the difference between the average energy of the total system and the average energy 
of the isolated leads can be written: 

&imp ~\~ ^imp— leads ~\~ &&leads (-^"-Q 

where: 

1 f d 

£im P = ^ Z) / 2^7 + ^ G ™ {UJ) (B " 2) 



Simp-leads ~ ^ / 7. — : ^0,nfc Cfcn (B"3) 



n,k 

(we made use of the fact that G < (u>y = —G < (u)), and: 



5S, 



leads 



The expression for G kn (u) can be derived rather easily in the present case of non- 



interacting leads 



49|: 



G£» = E G<» + £ y< ( W ) if , fcm (B-5) 

Using: 
and 

E -^°' nfe S'fc*^) -^o,fcm = ifL{u)T Lnm (uj) + i/ R (u;)r Rnm (u;) (B-7) 
one arrives at the following expression for £i mp -ieads'- 

Simp-leads =®J 7^ Tr {l( A l^) + A «M] G< H + i[h{u)T L {u) + f R {u)T R {u)) G\u)} 

(B-8) 

Using also the fact that the flux of particles coming in and out from the junction is exactly 
zero in steady states: 

J duo [J L (u) + J R {u)\ = (B-9) 

J dco Tr{[T L (u) + r R (u)) G<(u) - \f L (u)T L (u) + f R (u)T R {u)) [G a (cu) - G r (cu)]} = 

(B-10) 

one can ignore taking the real part of the r.h.s. of Eq. ( 1B-8I) : 

Simp-leads = J ^ Tr{[(A r L (uj) + A R (u)} G<(u) + i[f L (u)T L (u) + f R (u)T R (u)] G a (oo)} 

(B-ll) 

which can be further written as in Eq. (l3Tj) . 

Now let's focus on the expression for SSi ea d s - Similarly to Eq. (1B-5P one also has: 



= + H om G< k (u) + 9kk(") Hom G a nk (u) (B-12) 



Further use of: 

G nk{ U ) = E G nm( U ) H °> mk 9kk{ U ) + H 0,mk 9kk( U )^ (B-13) 
m m 

GIM = E G nm^) H 0,mk (B-14) 
m 

£fcfc iM M = (w) [9tk L(R) M - ^ fcL(fl) Ml , (B-15) 

£ + -) ^ <4H ^L(-) ^o, fe n = lim| £ { Jl^™r® iS) (B-16) 

feei(i?) ^ v ;v ; 

allows us to write the expression for 5£i ea d s as: 

= \f^l Tr{[S L (u)+S R (u)] G<(lu)-[S l (lu) f L (co)+S R (co) f R (u)} 



Tr{[F L (u;) f L {u) + F R {u) f R {u)\ G r (u)} + h.c] (B-17) 

with: 

S L (R) nm {u;) = lim / — — — — B-18 

s^o J 2tc [oj — e + to)(uj — e — i5) 

The function 5" (a;) has a singular part which however doesn't contribute to 5£\ ea d s - Indeed, 
writing: 

1 f 5 5 

ShruM = $tF L(R)rim (u;) + \im - J de (u + e) T L(R)nm {e) j- _ ^ - ^ ( ^_ e)2 + <52 , 

(B-20) 

the contribution to 5£i ea d s of the second term on the r.h.s. of Eq. (1B-20I) is proportional to: 
lim ~ / du u Tr{ [T L (u) +T R (u)] G < (w) - [f L (uj)r L (uj)+f R (uj)r R (uj)] [G a (u)-G r (u)}} = 

5^0 5 J 

(B-21) 

which vanishes by virtue of the fact that the integral multiplying | is proportional to the 
flux of energy coming in and out from the junction, which is exactly zero in steady states: 

J duo uj [J l {uj) + J r {oj)\ = (B-22) 

Thus, the expression for 5£i ea ds becomes: 

5£i eads = \J^l Tr{[^F L (uj)+^F R (uj)] G<(uj)-[^F l (uj) f L (uj)+^F R (uj) f R (u)] [G a (u)-G r (u)]} 



rr{[F £ («) / L (w) + F a (w) f R (u>)) G r (uj)} (B-23) 

Noting that the function i^rm (a>) is related in a simple way to the energy derivative of 
A£(m(ci>), one finally arrives at Eqs. 



APPENDIX C: AVERAGE NUMBER OF DISPLACED ELECTRONS IN THE 

LEADS 

In a manner similar to the one described in detail in Appendix B, one can obtain an 
expression for the average number of electrons displaced in the two leads: 

5N leads = E/^ " 9kkH\ (C-1) 

with the final expression reading: 

SN leads = - fp-. Trii^AKu) + ^A r R (u)} G<(u) 
J 2ix% aw aio 

- i[ h(co) ^A^) + f a {u>) ^QfA^H] [G a (co) + (C-2) 
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FIG. 1: Schematic view of the Anderson impurity model system considered, (a) Tight binding 
model for the non-interacting system, (b) Definition of applied source-drain bias V and gate 
voltage Vg- 




FIG. 2: Results for the self-consistent GW approximation at zero applied source-drain bias and 
ksT = 0. (a) Impurity occuption number as a function of gate voltage, (b) Spectral function for 
three different values of the gate voltage Vg- Using U = 4.78 and T = 0.05. 




FIG. 3: False color plots of junction properties calculated in the self-consistent GW approximation 
as a function of the applied source-drain bias V and gate voltage Vq at fceT = 0. (a) Current, (b) 
Average impurity occupation number. Using U = 4.78 and T = 0.05. 
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FIG. 4: Results for a non-magnetic solution throughout the gate bias range in the GW approx- 
imation at zero source-drain bias and ksT = 0.01. (a) Impurity occuption number as a function 
of gate voltage, (b) Spectral function for two values of gate voltage, the symmetric case and an 
asymmetric case. Using U = 4.78 and T = 0.05. 
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FIG. 5: (a) Current as a function of the applied bias for gate voltage fixed to the symmetric 
case and k^T = 0.01. Curves labeled magnetic correspond to a bias sweep from V=0 to V=8. 
Curves labeled non-magnetic correspond to a reverse bias sweep from V=8 to V=0. Results for 
the Hartree-Fock and GW approximations are compared, (b) Corresponding spectral functions for 
applied source-drain bias V = 2. (c) Comparison of spectral functions for the non-magnetic solution 
in the GW approximation at three different applied source-drain bias values. Using U = 4.78 and 

r = cos. 




FIG. 6: Real and imaginary parts of the retarded self-energy in the GW approximation for the 
non-magnetic solution at half-filling, applied source-drain bias V = 2 and k^T = 0.01. Using 
U = 4.78 and T = 0.05. 
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FIG. 7: Illustration of steps in the iterative solution to arrive at the final non-magnetic solution 
in the GW approximation. Gate voltage fixed to the symmetric (half-filling) case, applied source- 
drain bias V = 2 and ksT = 0.01. (a) Spectral function for the non-magnetic Hartree-Fock and 
GqWq solutions, (b) Real and imaginary parts of the retarded non-magnetic GoWo self-energy. 
Using U = 4.78 and T = 0.05. 




FIG. 8: Change in the average energy of the system as a function of applied source-drain bias for 
gate voltage fixed to the symmetric case (half-filling) and fc^T = 0.01. Curves labeled magnetic 
correspond to a bias sweep from V=0 to V=8. Curves labeled non-magnetic correspond to a 
reverse bias sweep from V=8 to V=0. Results for the Hartree-Fock and GW approximations are 
compared. Using U = 4.78 and T = 0.05. 



